A Monte Carlo investigation of the critical behavior of Stavskaya's probabilistic cellular automaton 
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Stavskaya's model is a one-dimensional probabilistic cellular automaton (PC A) introduced in the end of the 
1960's as an example of a model displaying a nonequilibrium phase transition. Although its absorbing state 
phase transition is well understood nowadays, the model never received a full numerical treatment to investigate 
its critical behavior. In this brief report we characterize the critical behavior of Stavskaya's PCA by means of 
Monte Carlo simulations and finite-size scaling analysis. The critical exponents of the model are calculated 
and indicate that its phase transition belongs to the directed percolation universality class of critical behavior, 
as it would be expected on the basis of the directed percolation conjecture. We also explicitly establish the 
relationship of the model with the Domany-Kinzel PCA on its directed site percolation line, a connection that 
seems to have gone unnoticed in the literature so far. 
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I. INTRODUCTION 

Stavskaya's model is a one-dimensional probabilistic cel- 
lular automaton (PCA) proposed in the end of the 1960's by 
the Russian school of Markov processes as an example of an 
interacting particle system presenting a nonequilibrium phase 
transition [1—3]. The model is related with the directed site 
percolation (DP) process, of which it can be viewed as a one- 
sided version, as well as with the Domany-Kinzel PCA in one 
of its manifolds (cf. Sec. II) [4]. Roughly speaking, the phase 
transition in Stavskaya's model follows from its attractiveness 
(its tendency for forming clusters) and the existence of an ab- 
sorbing state, and is well understood nowadays. 

However, while many rigorous results exist for this model 
[1-3, 5-11], it has never received a full numerical treatment 
to estimate its critical point and critical exponents. In this 
brief report we proceed to such an investigation of Stavskaya's 
model by Monte Carlo simulations and finite-size scaling 
techniques. Besides closing a gap in the characterization of 
the model, our results add another bit of evidence in favor of 
the DP conjecture, according to which phase transitions into 
an absorbing state in short-ranged single component systems 
in the absence of conserved quantities all belong to the same 
universality class of critical behavior [12, 13]. Remind that 
although the DP conjecture is grounded on solid theoretical 
arguments and has been verified in a host of model systems, it 
(/) could not be proved rigorously yet, and (ii) has only very 
thin experimental evidence [14, 15], so that it continues to rely 
on model systems to sustain itself. 

II. STAVSKAYA'S MODEL 

Stavskaya's model is a two-state PCA defined on a one- 
dimensional periodic lattice of L cells specified by the con- 
figuration 77(f) = (r?i(t), rj 2 {t), r) L (t)) € {0, 1} L and 
evolving in discrete time t E N according to the following 
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very simple rule: with probability e £ [0, 1], rji(t + 1) = 1, 
otherwise r)i(t + 1) = r)i-i(t) ■ T]i(t). 

Clearly, 1 = (1, 1, . . . , 1) is an absorbing state of the 
model. It can be proven that there exists a critical e* such that 
for e > e* the only invariant measure of Stavskaya's PCA is 
Si, the measure concentrated in 1, and that for e < e* the 
invariant measures are translation-invariant convex combina- 
tions of the form a>fj, e + (1 — a)Si, with < a < 1 and 
jU e the measure that puts mass on configurations with density 
< < 1 [2, 3, 6]. Early bounds on the critical point es- 

timate 0.09 < e* < 0.323 [3, 6]. The upper bound was even- 
tually confirmed, but not improved, by different techniques 
[8], while the lower bound never received a reassessment; it 
should be remarked that lower bounds on critical values of 
interacting particle systems are notoriously difficult to obtain. 

Stavskaya's PCA is related with the Domany-Kinzel (DK) 
PCA [4] by taking the complementary (negated) variables 
fji = 1 — rji- It can then be seen that Stavskaya's PCA corre- 
sponds to the DK PCA on the line p\ = p 2 = 1 — e, i.e., 
over the directed site percolation (site DP) line of the DK 
PCA parameter space. Notice, however, that the dynamics 
in the DK PCA is defined for each of its two sublattices in 
the time direction (even and odd time steps), while the dy- 
namics in Stavskaya's PCA is direct. On the site DP line, the 
DK PCA displays an inactive-active phase transition at the 
critical point p* = p\ = 0.705 489(4) [17], corresponding 
to e* = 0.294 511(4), within the rigorous bounds mentioned 
before (the numbers between parentheses indicate the uncer- 
tainty in the last digit or digits of the data). Curiously, the re- 
lationship between Stavskaya's PCA and the DK PCA seems 
to have gone unnoticed in previous investigations [6-11], al- 
though a coupling scheme with an "independent oriented per- 
colation" process equivalent with site DP was used in [8]. It 
is worth mentioning that Stavskaya's model, together with an- 
other PCA introduced by the same epoch, Vasil'ev's model 
[6, 18] — which corresponds to the p 2 = line in the DK PCA 
or, equivalently, to a probabilistic version of CA rule 18 in 
Wolfram's classification scheme [19] — predates the DK PCA 
and related models by almost two decades, but did not receive 
much attention, not even when CA and PCA reentered the 
mainstream scientific agenda in the 1980's. 
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FIG. 1. Stationary density of inactive cells in an automaton of L = 
8000 cells averaged over 1000 samples. The inset shows the curve 
close to the critical point e* = 0.29450(5) (value obtained from 
time-dependent simulations). 



III. THE CRITICAL BEHAVIOR 

Our Monte Carlo (MC) simulations of Stavskaya's model 
ran as follows. For a given s, the PC A is initialized with each 
rji(0) = 1, 1 ^ i ^ L, drawn independently with probability 
1/2. Stationary state quantities, e.g. the density of active cells 
Pl = L^ 1 Vi< are men sampled after the system is relaxed 
through 100L MC steps (MCS), with one MCS equivalent to 
a synchronous update of the states of all L cells of the au- 
tomaton. This amount of relaxation proved enough for our 
purposes. Moreover, except for the data shown in Figure 1, 
our results were obtained from time-dependent simulations, 
so that estimates on the stationary state did not concern us 
much. We refer the reader to [16] for a nice exposition of the 
time-dependent techniques employed in what follows. 

The critical behavior of the model can be determined by 
assuming the scaling relation 



1 - PL (t- A) ~ t-P/"ii $(A t^"w , t v ^' v w /L) (1) 

close to the critical point e* L , with A = e — e* L ^ 0. We 
do not put a subscript 'U on A or the critical exponents to 
lighten the notation. For a very large system, relation (1) be- 
comes 1 - pL(t;A) ~ t~^l v \\ $(Ai 1 /' y ii), with $(x < 1) ~ 



constant and $(x ^ 1) 



rP 



The investigation of the time- 



dependent profiles pi,(t;A) then allow for the simultaneous 
determination of e* L and 5 = ft/vn, and judicious perusal of 
(1) and derived relations furnish the other exponents. 

Figure 1 displays the density profile p^ for an automaton 
of L = 8000 cells in the stationary state. We actually plot the 
density of inactive cells 1 — pl instead, because it is this quan- 
tity that enters the scaling relation (1). The steep transition 
about e* L anticipates a small value for the exponent (3. To esti- 
mate e* more precisely, we plot 1 — close to s ~ 0.294 
for some large L. On the critical point, 1 — pl(^) ~ t~ s and 
we can estimate 6 by plotting log b [(l — pi(t/&))/(l — PL(t))] 
against 1/t for some small b. Our data for L = 20000 and 
b = 10 appear in Figure 2. From these data we could extract 
the estimates e* = 0.29450(5) and S = 0.155(5). Similar 
estimates using 4000 ^ L ^ 16000 confirm these values. 
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FIG. 2. (a) Logarithmic plot of 1 - p L (t) with L = 20000 cells. 
(b) Instantaneous values of S obtained from the curves in (a). In 
both graphs we have, from the uppermost curve downwards, e = 
0.2943, 0.2944, 0.2945, 0.2946, and 0.2947. From these curves we 
estimated e* = 0.29450(5) and 5 = 0.155(5). The dashed line in 
panel (6) indicates the best value available for <5dp- 



Notice that this estimate of e* completely agrees with the crit- 
ical point p\ = pl = 0.705 489(4) = 1 - e* found for the site 
DP transition in the DK PC A [17]. 

The exponent un can be obtained by plotting t s (l — Pb{t)) 
versus iA^ii and tuning v\\ to achieve data collapse with dif- 
ferent A. The collapsed curves shown in Figure 3 were ob- 
tained with a combination of central values e* — 0.29451, 
5 = 0.157, and v\\ = 1.73. We found it hard to discern values 
of vu by less than ±0.02. Otherwise, we found the data col- 
lapse very sensitive to the choice of e* ; in fact, it could have 
been used to locate e* within quite tight bounds. Combining 
S and v\\ furnishes /3 = 5v\\ = 0.27(1) (or /? = 0.268(10)). 

The third independent exponent can be obtained by plotting 
t s (l — pL{t)) versus t/ L z for different L and tuning z until 
data collapse for some z. Since z = v\\/v± by definition, 
this procedure also gives v±_ once Mi is known. The finite-size 
curves appear in Figure 4. We found z = 1.6(1). These three 
exponents, 5, vn, and z, suffice to determine the universality 
class of critical behavior of the model, the other exponents 
following from well known hyperscaling relations [16]. 

The best values available for 5, u», (3, and z for the 
DP process on the square lattice are <5rjp = 0.159 464(6), 
^IIdp = 1.733 847(6), f3 DP = 0.276 486(8), and z DP = 
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FIG. 3. Data collapse of the scaled time-dependent density profiles 
for e — e* = ±0.0001, ±0.0002. The upper (lower) branches corre- 
spond to e < e* (e > e*). The best data collapse was obtained with 
the central values e* = 0.29451, 8 = 0.157, and v« = 1.73. 




FIG. 4. Finite-size data collapse of the scaled time-dependent density 
profiles on the critical point e* = 0.29450 for 4000 < L < 20000. 
Best data collapse was obtained with S = 0.157 (the same value as in 
Fig. 3) and z — 1.6. Notice that the data is spread over > 3 decades 
and that the range in the y-axis is reasonably tight. 



1.580 745(10) [20]. Thus, within the error bars our estimates 
for these exponents put the critical behavior of Stavskaya's 
model phase transition in the DP universality class, as it would 
be expected on the basis of the DP conjecture. 



IV. SUMMARY AND CONCLUSIONS 

We estimated the critical point of Stavskaya's model at 
e* = 0.29450(5) and found that the model belongs to the 
DP universality class of critical behavior. The value of 1 — e* 
is in excellent agreement with the critical point p\ = p* 2 = 
0.705 489(4) for the site DP transition in the DK PCA [17]. 

The estimates of the critical point as well as of the criti- 
cal exponents could be improved by larger simulations, but 
we believe that this would be superfluous, since both the lo- 
cation of the critical point of the model within better bounds 
than those provided by rigorous and mean-field analyses and 
the determination of its universality class of critical behavior 
could be established within the computational efforts reported 
here, namely, a few thousand hours of CPU time on Intel i7- 
860 processors running GCC/Linux at 2.8 GHz. 

Taken together, our numerical results for Stavskaya's PCA 
and the rigorous results existent on its relationship with the 
general theory of cellular automata and percolation processes 
provide a reasonably complete characterization of the model. 
The establishment of its relationship with the DK PCA on 
its site DP line allows the translation of results between the 
two models, with potential benefits for future developments 
involving either model per se or as approximations (e.g., on 
coupling schemes) to models of greater complexity. 
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